#include "ImgQuaEv.h"

using namespace std;

bool ImageQuaEv::IsValidInts(int imgarr[])
{
	if (imgarr[1] > imgarr[3])
	{
		if (imgarr[0] > imgarr[3] || imgarr[4] > imgarr[0])
		{
			return false;
		}
	}
	else
	{
		if (imgarr[1] < imgarr[4] || imgarr[4] < imgarr[0])
		{
			return false;
		}
	}
	return true;
}

namespace ImageQuaEv
{	
	float modeInts(const cv::Mat&Img512, int&modeInts) //return MeanInts
	{
		uint64 sum = 0;
		std::vector<ushort>histDis(65536, 0);
		ushort value_pixel = 0;
		int R = Img512.rows;
		int C = Img512.cols;
		if (Img512.type() == CV_8UC1)
		for (int i = 0; i < R; i++)
		{
			const uchar *p = Img512.ptr<uchar>(i);
			for (int j = 0; j < C; j++)
			{
				value_pixel = *(p + j);
				sum += value_pixel;
				histDis[value_pixel]++;
			}
		}
		else
		for (int i = 0; i < R; i++)
		{
			const ushort *p = Img512.ptr<ushort>(i);
			for (int j = 0; j < C; j++)
			{
				value_pixel = *(p + j);
				sum += value_pixel;
				histDis[value_pixel]++;
			}
		}
		float mInts = 1.0f*sum / R / C;
		//////////////
		modeInts = 0;
		int maxCount = 0;
		for (int i = 0; i < 65536; i++)
		{
			if (histDis[i]>maxCount)
			{
				maxCount = histDis[i];
				modeInts = i;
			}
		}
		return mInts;
	}

	template<typename T>
	bool isBasicPeak(const int x, const int y, cv::Mat& img, const int threshold)
	{
		int center = img.at<T>(x, y);
		int pixel3_3[9] = { img.at<T>(x - 1, y - 1),
			img.at<T>(x - 1, y),
			img.at<T>(x - 1, y + 1),
			img.at<T>(x, y - 1),
			center,
			img.at<T>(x, y + 1),
			img.at<T>(x + 1, y),
			img.at<T>(x + 1, y - 1),
			img.at<T>(x + 1, y + 1)
		};
		int falseNum = 0;
		for (int i = 0; i < 4; i++)
		{
			if (center < pixel3_3[i]) return false;
			if (center == pixel3_3[i]) falseNum += 1;
			if (threshold > pixel3_3[i]) return false;
		}
		for (int i = 5; i < 9; i++)
		{
			if (center < pixel3_3[i]) return false;
			if (center == pixel3_3[i]) falseNum += 1;
			if (threshold > pixel3_3[i]) return false;
		}
		if (falseNum > 2) return false;
		return true;
	}

	bool basicPeak(const int x, const int y, cv::Mat& img, const int threshold)
	{
		return isBasicPeak<unsigned short>(x, y, img, threshold);
	}
	bool basicPeak8bit(const int x, const int y, cv::Mat& img, const int threshold)
	{
		return isBasicPeak<uchar>(x, y, img, threshold);
	}

	template<typename T>
	bool isGoodPeak(int x, int y, cv::Mat& img)
	{
		int count = 0;
		T centerVal = img.at<T>(x, y);
		if (img.at<T>(x - 1, y - 1) < img.at<T>(x - 2, y - 2)) count++;
		if (img.at<T>(x + 1, y + 1) < img.at<T>(x + 2, y + 2)) count++;
		if (img.at<T>(x + 1, y - 1) < img.at<T>(x + 2, y - 2)) count++;
		if (img.at<T>(x - 1, y + 1) < img.at<T>(x - 2, y + 2)) count++;
		if (img.at<T>(x - 1, y) < img.at<T>(x - 2, y)) count++;
		if (img.at<T>(x + 1, y) < img.at<T>(x + 2, y)) count++;
		if (img.at<T>(x, y - 1) < img.at<T>(x, y - 2)) count++;
		if (img.at<T>(x, y + 1) < img.at<T>(x, y + 2)) count++;

		int a[5] = { 0 };
		a[0] = img.at<T>(x - 2, y);
		a[1] = img.at<T>(x - 1, y);
		a[2] = centerVal;
		a[3] = img.at<T>(x + 1, y);
		a[4] = img.at<T>(x + 2, y);
		if (!IsValidInts(a)) count++;
		a[0] = img.at<T>(x, y - 2);
		a[1] = img.at<T>(x, y - 1);
		a[3] = img.at<T>(x, y + 1);
		a[4] = img.at<T>(x, y + 2);
		if (!IsValidInts(a)) count++;
		return count < 3;
	}

	bool goodPeakJudge8bit(int x, int y, cv::Mat& img, const int threshold)
	{
		if (!basicPeak8bit(x, y, img, threshold))
		{
			return false;
		}

		return isGoodPeak<unsigned char>(x, y, img);
	}
	bool goodPeakJudge(int x, int y, cv::Mat& img, const int threshold)
	{
		if (!basicPeak(x, y, img, threshold))
		{
			return false;
		}
		return isGoodPeak<unsigned short>(x, y, img);
	}

	int getMultxy(int x, int y, const cv::Mat& img, float threshold, bool is16Bit)
	{
		auto getInts = [&](int xL, int yL) {
			unsigned short Ints = is16Bit ? img.at<unsigned short>(xL, yL) : img.at<unsigned char>(xL, yL);
			return Ints;
			};
		unsigned short CenterInts = getInts(x, y);// isBit16 ? img.at<unsigned short>(x, y) : img.at<unsigned char>(x, y);
		unsigned int sum = 0;
		int c = img.cols;
		int r = img.rows;
		int xLen = 1;
		while (x + xLen < r)
		{
			unsigned short Ints = getInts(x + xLen, y);
			if (Ints < threshold)
				break;
			if (Ints / CenterInts > 1)
				sum++;
			if (sum > 3)
				return 0;
			xLen++;
		}

		int p = 0;
		p = 1; sum = 0;
		while (x - p > 0)
		{
			unsigned short Ints = getInts(x - p, y);
			if (Ints < threshold)
				break;
			if (Ints / CenterInts > 1)
				sum++;
			if (sum > 3)
				return 0;
			p++;
		}
		//int rateX = p>xLen ? p / xLen : xLen / p;
		if ((p > xLen ? p / xLen : xLen / p) > 2)
			return 0;
		xLen += p - 1;

		sum = 0;
		int yLen = 1;
		while (y + yLen < c)
		{
			unsigned short Ints = getInts(x, y + yLen);
			if (Ints < threshold)
				break;
			if (Ints / CenterInts > 1)
				sum++;
			if (sum > 3)
				return 0;
			yLen++;
		}
		p = 1; sum = 0;
		while (y - p > 0)
		{
			unsigned short Ints = getInts(x, y - p);
			if (Ints < threshold)
				break;
			if (Ints / CenterInts > 1)
				sum++;
			if (sum > 3)
				return 0;
			p++;
		}
		if ((p > yLen ? p / yLen : yLen / p) > 2)
			return 0;
		yLen += p - 1;
		return xLen * yLen;
	}
	int getXmultiY(int x, int y, const cv::Mat& img, float threshold)
	{
		return getMultxy(x, y, img, threshold, true);
	}
	int getXmultiY8bit(int x, int y, const cv::Mat& img, float threshold)
	{
		return getMultxy(x, y, img, threshold, false);
	}

	float getScore(int x, int y, const cv::Mat &img, int s = 1)
	{
		int EightPixelInts = img.at<unsigned short>(x - s, y) +
			img.at<unsigned short>(x - s, y - s) + img.at<unsigned short>(x - s, y + s) +
			img.at<unsigned short>(x + s, y) + img.at<unsigned short>(x + s, y - s) +
			img.at<unsigned short>(x + s, y + s) + img.at<unsigned short>(x, y - s) +
			img.at<unsigned short>(x, y + s);
		int centerPixelInts = img.at<unsigned short>(x, y);
		float meanPeakInts = 1.0*(EightPixelInts + centerPixelInts) / 9;
		float score = 1.0*(8 * centerPixelInts - EightPixelInts) / meanPeakInts;

		return score;
	}

	float divideBack(cv::Mat&image)//not used
	{
		cv::Mat bg;
		cv::Mat element = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(30, 30));
		morphologyEx(image, bg, cv::MORPH_OPEN, element);
		cv::subtract(image, bg, image);
		float m = cv::mean(image)[0];
		return m;
	}

	template<typename T>
	void isValidPeak(const cv::Mat& img, std::vector<cv::Point2i>& peak, bool is16Bit)
	{
		peak.clear();
		std::vector<cv::Point2i> peakTemp;
		peakTemp.reserve(10000);
		std::vector<unsigned short> intRecord;
		intRecord.reserve(10000);
		int c = img.cols;
		int r = img.rows;
		//////////////		
		cv::Mat tempImg = img.clone();
		//float mInts = divideBack(tempImg);
		int modeIntsVal = 0;
		float mInts = modeInts(tempImg, modeIntsVal);
		if (mInts < 1)	return;
		int p12 = mInts * 2; //2.3  2
		int p11 = mInts * 1.4; //2 1.4
		int p20 = mInts * 4;
		for (int x = 5; x < r - 5; x++)
			for (int y = 5; y < c - 5; y++)
			{
				unsigned short Ints = tempImg.at<T>(x, y);
				bool peakFlag = false;
				if (Ints > p12) //1.4
				{
					if (Ints > p20)
					{
						peakFlag = is16Bit ? basicPeak(x, y, tempImg, p11) : basicPeak8bit(x, y, tempImg, p11);
					}
					else
					{
						peakFlag = is16Bit ? goodPeakJudge(x, y, tempImg, p11) : goodPeakJudge8bit(x, y, tempImg, p11);
					}

					if (peakFlag) //find peak 1.1
					{
						int area = is16Bit ? getXmultiY(x, y, tempImg, p11 * 0.8) : getXmultiY8bit(x, y, tempImg, p11 * 0.8);
						if (area) //ret  trues
						{
							peakTemp.push_back(cv::Point2i(y, x));
							intRecord.push_back((int)img.at<uchar>(x, y));
						}

					}
					//if (peakFlag) //find peak 1.1
					//{
					//	int area = is16Bit ? getXmultiY(x, y, tempImg, p11 * 0.8) : getXmultiY8bit(x, y, tempImg, p11 * 0.8);
					//	if (area)
					//	{
					//		peak.push_back(cv::Point2i(y, x));
					//	}
					//	intRecord.push_back(img.at<T>(x, y));
					//}
				}
			}

		std::vector< unsigned short >::iterator p98 = intRecord.begin() + intRecord.size() * 98 / 100; //p90
		std::vector< unsigned short >::iterator p80 = intRecord.begin() + intRecord.size() * 80 / 100; //p90
		std::nth_element(intRecord.begin(), p98, intRecord.end());
		std::nth_element(intRecord.begin(), p80, intRecord.end());

		int largeThresh = *p98;
		int smallThresh = *p80;
		cv::Mat b_m, b_st;
		cv::meanStdDev(tempImg, b_m, b_st);
		smallThresh = 30 / (2 * b_st.at<double>(0, 0) + 0.1);
		for (int pointNum = 0; pointNum < intRecord.size(); pointNum++) {
			if (intRecord[pointNum] > smallThresh && intRecord[pointNum] <= largeThresh) {
				peak.push_back(peakTemp[pointNum]);
			}
		}

		////过曝临时方案-chenwei
		//for (int pointNum = 0; pointNum < intRecord.size(); pointNum++) {
		//	peak.push_back(peakTemp[pointNum]);
		//}
	}

	void SerachGoodPeak16bit(const cv::Mat &img, std::vector<cv::Point2i> &peak)//score peak locate
	{
		isValidPeak<unsigned int>(img, peak, true);
	}
	void SerachGoodPeak8bit(const cv::Mat &img, std::vector<cv::Point2i> &peak)//score peak locate
	{
		isValidPeak<uchar>(img, peak, false);
	}

	template<typename T>
	void getScoreSlop(int x, int y, const cv::Mat& img, std::vector<float>& intentRecord)
	{
		int ids = 0;
		for (int i = y - 5; i < y + 6; i++) {
			intentRecord[ids] += img.at<T>(x, i);
			ids += 1;
		}
	}

	
	int CalcMatQua(const cv::Mat &img, double &midScore, double &intensity, float intenrecord_ret[11] = nullptr) //return peakNumber
	{		
		std::vector<float> intensity_record(11, 0.0);
		midScore = 0;

		std::vector<cv::Point2i> peak;
		CV_16UC1 == img.type() ? SerachGoodPeak16bit(img, peak) : SerachGoodPeak8bit(img, peak);
		if (peak.size() < 9)
		{
			return 0;
		}

		for (uint i = 0; i < peak.size(); i++)
		{
			cv::Point2i& p = peak[i];
			ImageQuaEv::getScoreSlop<uchar>(p.y, p.x, img, intensity_record);
		}
        
        float min_intent = 153000;
		std::for_each(intensity_record.begin(), intensity_record.end(), [&](float& inten_v) {
			inten_v /= peak.size();
			min_intent = (min_intent > inten_v ? inten_v : min_intent);
		});

        intensity = intensity_record[5];
        float meanIntent = (intensity_record[5] + intensity_record[3] + intensity_record[7]) / 3;
        float quavalue = (intensity_record[5] * 2 - intensity_record[3] - intensity_record[7]) / meanIntent;

  //      float sum = 0;
		//std::for_each(intensity_record.begin(), intensity_record.end(), [&](float& inten_v) {inten_v -= min_intent; sum += inten_v; });
  //      sum /= intensity_record.size();

		float sum = 0;
		for (auto i = 0; i < intensity_record.size(); i++) {
			intensity_record[i] -= min_intent;
			sum += intensity_record[i];
			if (intenrecord_ret)
			{
				intenrecord_ret[i] = intensity_record[i];
			}
		}
		sum /= 11;


        float stds = 0;
		std::for_each(intensity_record.begin(), intensity_record.end(), [&](float& inten_v) { stds += std::pow(inten_v - sum, 2); });
		stds /= intensity_record.size();

        midScore = (stds / sum + 40 * quavalue) / 4;
		return (int)peak.size();		
	}

	int GetImgQua(const cv::Mat &img, double &midScore, double &intensity, int nRoiSize)
	{
		int R = img.rows;
		int C = img.cols;
		int midSize = nRoiSize;
		int startPosiY = R / 2;
		int startPosiX = C / 2;
		if (R > midSize) R = midSize;
		if (C > midSize) C = midSize;

		cv::Rect cutImgSize(startPosiX - C / 2, startPosiY - R / 2, C, R);
		cv::Mat cutImg = img(cutImgSize);
       int peaksize = CalcMatQua(cutImg, midScore, intensity);
		midScore /= 2;
		return peaksize;
	}

	//std::vector<float>& intenvec, int nRoiSize = 256
	int GetImgDetialQua(const cv::Mat& img, double& midScore, double& intensity, std::vector<float>& intenvec, int nRoiSize)
	{
		int R = img.rows;
		int C = img.cols;
		int midSize = nRoiSize;
		int startPosiY = R / 2;
		int startPosiX = C / 2;
		if (R > midSize) R = midSize;
		if (C > midSize) C = midSize;

		cv::Rect cutImgSize(startPosiX - C / 2, startPosiY - R / 2, C, R);
		cv::Mat cutImg = img(cutImgSize);
		float intensityarr[11];
		int peaksize = CalcMatQua(cutImg, midScore, intensity, intensityarr);
		midScore /= 2;
		for (auto i = 0; i < 11; i++)
		{
			intenvec.push_back(intensityarr[i]);
		}
		return peaksize;
	}

	int midScoreSelf(const cv::Mat& img, double& midScore, double& intensity)
	{
		int peaksize = CalcMatQua(img, midScore, intensity);
		midScore /= 2;
		return peaksize;
	}


	std::vector<std::vector<int>>  ScoreHeatMapSelf(const cv::Mat& totalImg, std::vector<std::vector<double>>& totalScore, std::vector<std::vector<double>>& totalInts)
	{
		int gridsize_x = 8;
		int gridsize_y = 4;
		int sw = int(totalImg.cols / gridsize_x);
		int sh = int(totalImg.rows / gridsize_y);
		int halfH = std::min(std::min(sw / 2, sh / 2), 512 / 2);
		int halfW = std::min(std::min(sw / 2, sh / 2), 512 / 2);
		double midScore, intensity;
		int sth, stw;
		cv::Mat subimg;
		std::vector<std::vector<int>> totalPeaksize;
		totalScore.resize(gridsize_y);
		totalInts.resize(gridsize_y);
		totalPeaksize.resize(gridsize_y);
		for (size_t iy = 0; iy < gridsize_y; ++iy) {
			totalScore[iy].resize(gridsize_x);
			totalInts[iy].resize(gridsize_x);
			totalPeaksize[iy].resize(gridsize_x);
			for (size_t ix = 0; ix < gridsize_x; ++ix) {
				if (halfH < 10 || halfW < 10) {
					totalScore[iy][ix] = 0.f;
					totalInts[iy][ix] = 0.f;
					totalPeaksize[iy][ix] = 1;
				}
				else {
					sth = sh*iy + sh / 2 - halfH;
					stw = sw*ix + sw / 2 - halfW;
					subimg = totalImg(cv::Rect(stw, sth, 2 * halfW, 2 * halfH));
					int peaksize = CalcMatQua(subimg, totalScore[iy][ix], totalInts[iy][ix]);
					totalScore[iy][ix] /= 2;
					totalPeaksize[iy][ix] = peaksize;
				}

			}
		}

		return totalPeaksize;
		
	}

	int midScoreSelfKK(const cv::Mat& img, double& midScore, double& intensity, bool &overDensity, bool &lessPoints, bool &overExposure)
	{
		int peaksize = CalcMatQuaKK(img, midScore, intensity, overDensity, lessPoints, overExposure);
		midScore /= 2;
		return peaksize;
	}

	void getFimg(std::vector<cv::Point2i>& resultPos, cv::Mat& img1, int __ImageHeight, int __ImageWidth) {

		int x, y;
		int spotwin = 3;
		const int fwhm = spotwin / 2;
		int csz = spotwin*sizeof(float);
		int k1, k2;
		int k1b, k1e, k2b;
		float rgmt[3][3] = { { 1.0f, 1.0f, 1.0f },{ 1.0f, 1.0f, 1.0f },{ 1.0f, 1.0f, 1.0f } };
		cv::Mat img0 = cv::Mat::zeros(__ImageHeight, __ImageWidth, CV_32F);
		for (size_t i = 0; i < resultPos.size(); ++i)
		{
			y = resultPos[i].y;
			x = resultPos[i].x;
			if (y > 6 && y < __ImageHeight - 6 && x > 6 && x < __ImageWidth - 6)
			{
				k1b = y - fwhm, k1e = y + fwhm, k2b = x - fwhm;
				for (k1 = k1b; k1 <= k1e; ++k1)
				{
					memcpy(img0.ptr<float>(k1) + k2b, rgmt[k1 - k1b], csz);
				}
			}
			img1.at<float>(y, x) = img1.at<float>(y, x) + 20;

		}

	}

	int CalcMatQuaKK(const cv::Mat &img, double &midScore, double &intensity, bool &overDensity, bool &lessPoints, bool &overExposure)
	{
		std::vector<float> intensity_record(11, 0.0);
		midScore = 0;
		std::vector<cv::Point2i> peak;

		std::vector<unsigned short> intRecord;
		intRecord.reserve(10000);

		//CV_16UC1 == img.type() ? SerachGoodPeak16bit(img, peak) : SerachGoodPeak8bit(img, peak);
		cv::Mat imag, noise;
		img.convertTo(imag, CV_16S);
		//
		std::vector<int> histData;
		int v1;
		histData.assign(65536, 0);
		float fillDensity = -99.0;
		fillDensity = getThresholdKK(imag, histData, 0, 0, 1);
		//
		int _ImageHeight = imag.rows, _ImageWidth = imag.cols;
		int sth;
		SpotFilterKK(imag, img, noise);
		OtsuThresholdKK(imag, sth, noise);
		fastLaunchKK(noise, imag, img, peak, sth);

		if (peak.size() < 9)
		{
			return 0;
		}

		for (uint i = 0; i < peak.size(); i++)
		{
			cv::Point2i& p = peak[i];
			ImageQuaEv::getScoreSlop<uchar>(p.y, p.x, img, intensity_record);
		}
		float min_intent = 153000;
		std::for_each(intensity_record.begin(), intensity_record.end(), [&](float& inten_v) {
			inten_v /= peak.size();
			min_intent = (min_intent > inten_v ? inten_v : min_intent);
		});
		// overDensity, lessPoints
		if (min_intent >= 240) overExposure = true;
		if (fillDensity > 0 && fillDensity < 0.04) overDensity = true;
		if (peak.size() < 20 && !overDensity) lessPoints = true;
		//

		intensity = intensity_record[5];
		float meanIntent = (intensity_record[5] + intensity_record[3] + intensity_record[7]) / 3;
		float quavalue = (intensity_record[5] * 2 - intensity_record[3] - intensity_record[7]) / meanIntent;

		//      float sum = 0;
		//std::for_each(intensity_record.begin(), intensity_record.end(), [&](float& inten_v) {inten_v -= min_intent; sum += inten_v; });
		//      sum /= intensity_record.size();

		float sum = 0;
		for (auto i = 0; i < intensity_record.size(); i++) {
			intensity_record[i] -= min_intent;
			sum += intensity_record[i];
			//if (intenrecord_ret)
			//{
			//	intenrecord_ret[i] = intensity_record[i];
			//}
		}
		sum /= 11;


		float stds = 0;
		std::for_each(intensity_record.begin(), intensity_record.end(), [&](float& inten_v) { stds += std::pow(inten_v - sum, 2); });
		stds /= intensity_record.size();

		midScore = (stds / sum + 40 * quavalue) / 4;
		return (int)peak.size();
	}

	void SpotFilterKK(cv::Mat& imag, const cv::Mat& rawim, cv::Mat& inoise) {
		cv::Mat bg, inoise_;
		float bgv, thv;
		cv::Mat m_kernel = (cv::Mat_<float>(3, 3) << -0.75, -0.5, -0.75,
			-0.50, 5.25, -0.50,
			-0.75, -0.5, -0.75);
		cv::Mat element = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(17, 17));
		morphologyEx(rawim, bg, cv::MORPH_OPEN, element);
		bg.setTo(cv::Scalar(1), bg == 0);
		bg.convertTo(bg, CV_16S);
		imag -= bg;
		bgv = cv::mean(bg).val[0];
		//thv = cv::mean(imag).val[0];
		cv::divide(imag, bg, inoise, 1, 5);
		int mp;
		auto thdv = binaryThresholdKK(inoise, imag, mp);
		inoise = inoise > thdv;
		inoise_ = imag > mp;
		cv::bitwise_and(inoise, inoise_, inoise);
		inoise.setTo(1, inoise);

		cv::resize(imag, imag, cv::Size(round(imag.cols * 1.25), round(imag.rows * 1.25)));
		imag.convertTo(imag, CV_32F);
		cv::GaussianBlur(imag, imag, cv::Size(3, 3), 0.9, 0.9);
		cv::filter2D(imag, imag, -1, m_kernel);

	}

	void OtsuThresholdKK(const cv::Mat& imgGrayt, int& sth, const cv::Mat& bimt) {

		cv::Mat imgGray = imgGrayt;
		cv::Mat bim = bimt;

		int i, j, sz;
		std::vector<int> histData;
		histData.assign(65536, 0);
		const int step = 256;
		sth = getThresholdKK(imgGray, histData, 1, step, 0) - step;

	}

	float binaryThresholdKK(const cv::Mat& imaget, const cv::Mat& imgt, int& mp) {

		cv::Mat image = imaget;
		cv::Mat img = imgt;
		std::vector<int> histData;

		double mv = cv::mean(image)[0];
		mv = std::max(mv, 0.01);
		histData.assign(65536, 0);
		int mst = 32;
		mp = getThresholdKK(img, histData, 0, mst, 0);
		int j;
		double sum = 0, sum_ = 0;
		int ssz = histData.size();
		int k = 0;
		for (j = 0; j < ssz; ++j) {
			sum += histData[j];
			if (j == mp) {
				sum_ = sum;
			}
		}
		// nead sum>0 and sum_>0;
		float rt = (0.6 + sum_ / sum);

		int i;
		float meanv = 0, stdv = 0;
		for (i = 0; i <= mp; ++i) {
			meanv += histData[i] * i;
		}
		meanv /= sum_;
		float meanvTmp = 0.f;
		for (i = 0; i <= mp; ++i) {
			meanvTmp = (i - meanv);
			stdv += histData[i] * meanvTmp * meanvTmp;
		}
		stdv = std::sqrt(stdv / sum_);

		auto temp = meanv + stdv * 3;
		if (mp < temp) mp = temp;
		mp -= 32;
		return mv*rt;
	}

	float getThresholdKK(cv::Mat& imgt, std::vector<int> &histData, int imageType, int offset, bool flagFillDensity) {

		int i, j, sz;
		cv::Mat img;
		int v1;
		img = imgt;
		//if (imageType)
		//	img = imgt + float(offset);
		//else
		//	img = imgt + short(offset);
		if (imageType) {
			for (size_t iy = 0; iy < img.rows; ++iy) {
				for (size_t ix = 0; ix < img.cols; ++ix) {
					v1 = img.at<float>(iy, ix);
					v1 += offset;
					if (v1 >= 0 && v1 < 65535) histData[v1]++;
				}
			}
		}
		else {
			for (size_t iy = 0; iy < img.rows; ++iy) {
				for (size_t ix = 0; ix < img.cols; ++ix) {
					//v1 = img.at<unsigned int>(iy, ix);
					v1 = img.at<short>(iy, ix);
					v1 += offset;
					if (v1 >= 0 && v1 < 65535) histData[v1]++;
				}
			}
		}
		int mx = -100;
		int mp = 0;
		const static int windows = 10;
		bool fg;
		sz = histData.size() - windows;
		int iu, id;
		int mn = (img.rows*img.cols) / histData.size() + 0.1;
		for (i = 0; i < windows; ++i) {
			if (mx < histData[i]) {
				mx = histData[i];
				mp = i;
			}
		}
		for (i = windows; i < sz; ++i)
		{
			if (histData[i] > mn)
			{
				iu = i + windows;
				id = i - windows;
				fg = true;
				for (j = id; j <= iu; ++j)
				{
					if (histData[i] < histData[j])
					{
						fg = false;
						break;
					}
				}
				if (fg)
				{
					mx = histData[i];
					mp = i;
					break;
				}
			}
		}
		//if (mx < mn) abort();
		if (mx < mn) {
			mx = histData[int(windows/2)];
			mp = int(windows/2);
		}
		//
		float fillDensity = (mx + 0.0001) / (img.cols*img.rows + 0.01);
		if (flagFillDensity) return fillDensity;
		const static float rt = 0.7;
		int mx_ = mx*rt;
		int mp_ = 0;
		sz = histData.size();
		for (i = mp; i < sz; ++i)
		{
			if (histData[i] < mx_)
			{
				mp_ = i;
				break;
			}
		}
		//if (mp_ < mp) abort();
		if (mp_ < mp) mp_ = mp+2;
		return mp + (mp_ - mp)*2.0;

	}


	void fastLaunchKK(cv::Mat& noise, cv::Mat& intensity, const cv::Mat img, std::vector<cv::Point2i>& resultPos, float sth) {

		cv::Mat imgmx, noise_;
		float itmp;
		float thdv = sth;
		int ib = 6, jb = 6;
		int ie = 0, je = 0;
		cv::Point2i pos;
		uchar num = std::floor(0.4 * 9.f) + 1;
		ie = intensity.rows - 6;
		je = intensity.cols - 6;
		cv::boxFilter(noise, noise_, -1, cv::Size(3, 3), cv::Point(-1, -1), false);
		cv::resize(noise_, noise_, cv::Size(round(noise.cols * 1.25), round(noise.rows * 1.25)), 0, 0, 1);
		cv::dilate(intensity, imgmx, cv::getStructuringElement(cv::MORPH_RECT, cv::Size(3, 3)));
		uchar* tptr1 = NULL;
		float* tptr2 = NULL;
		//
		resultPos.clear();
		std::vector<cv::Point2i> peakTemp;
		peakTemp.reserve(10000);
		std::vector<unsigned short> intRecord;
		intRecord.reserve(10000);
		//
		for (size_t i = ib; i < ie; ++i)
		{
			tptr1 = noise_.ptr<uchar>(i);
			tptr2 = intensity.ptr<float>(i);
			for (size_t j = jb; j < je; ++j)
			{
				itmp = tptr2[j];
				if (itmp > thdv && tptr1[j] > num)
				{
					if (imgmx.at<float>(i, j) != itmp) continue;
					pos.x = (int)(j/1.25 + 0.495);
					pos.y = (int)(i/1.25 + 0.495);
					peakTemp.push_back(pos);
					intRecord.push_back((int)img.at<uchar>(pos.y, pos.x));
				}
			}
		}
		std::vector< unsigned short >::iterator p98 = intRecord.begin() + intRecord.size() * 98 / 100;
		std::vector< unsigned short >::iterator p80 = intRecord.begin() + intRecord.size() * 80 / 100;
		std::nth_element(intRecord.begin(), p98, intRecord.end());
		std::nth_element(intRecord.begin(), p80, intRecord.end());
		int largeThresh = *p98;
		int smallThresh = *p80;
		for (int pointNum = 0; pointNum < intRecord.size(); pointNum++) {
			if (intRecord[pointNum] > smallThresh && intRecord[pointNum] <= largeThresh) {
				resultPos.push_back(peakTemp[pointNum]);
			}
		}
	}



}